HW5
Mike Cai
2024-05-30
1) Problem 14.3
#import data
data <- read_excel("Chapter14_exercises_data.xls", sheet = "Exercise 3, 4, 10")
sp500_ts <- ts(data$`Adj Close`, start = 2000, freq = 252)
plot(sp500_ts)
There is noticeable high volatility around the early 2000s and 2008,
corresponding to the dotcom bubble and the financial crisis,
respectively. These periods were marked by significant instability and
substantial fluctuations in stock prices. In contrast, there were also
calmer periods of lower volatility, such as the stretch from 2003 to
2007.
#get daily returns
returns <- diff(log(sp500_ts))
tsdisplay(returns, lag.max = 40)
returns_sq <- returns^2
tsdisplay(returns_sq, lag.max = 40)
The returns data seem to exhibit white noise, with the acf/pacf plots
not really showing decay or significant autocorrelation. The amplitude
of 0.05 is quite low as well. After squaring the data and looking at the
acf/pacf plots, there is decay and high persistence in the ACF, with
spikes cutting off at lag 14 in the PACF. I will try fitting an
ARCH(14).
model=ugarchspec(
variance.model = list(model = "sGARCH", garchOrder = c(14, 0)),
mean.model = list(armaOrder = c(2, 2), include.mean = TRUE),
distribution.model = "sstd")
#sanity check: explore the model parameters
model@model$pars## Level Fixed Include Estimate LB UB
## mu 0 0 1 1 NA NA
## ar1 0 0 1 1 NA NA
## ar2 0 0 1 1 NA NA
## ma1 0 0 1 1 NA NA
## ma2 0 0 1 1 NA NA
## arfima 0 0 0 0 NA NA
## archm 0 0 0 0 NA NA
## mxreg 0 0 0 0 NA NA
## omega 0 0 1 1 NA NA
## alpha1 0 0 1 1 NA NA
## alpha2 0 0 1 1 NA NA
## alpha3 0 0 1 1 NA NA
## alpha4 0 0 1 1 NA NA
## alpha5 0 0 1 1 NA NA
## alpha6 0 0 1 1 NA NA
## alpha7 0 0 1 1 NA NA
## alpha8 0 0 1 1 NA NA
## alpha9 0 0 1 1 NA NA
## alpha10 0 0 1 1 NA NA
## alpha11 0 0 1 1 NA NA
## alpha12 0 0 1 1 NA NA
## alpha13 0 0 1 1 NA NA
## alpha14 0 0 1 1 NA NA
## beta 0 0 0 0 NA NA
## gamma 0 0 0 0 NA NA
## eta1 0 0 0 0 NA NA
## eta2 0 0 0 0 NA NA
## delta 0 0 0 0 NA NA
## lambda 0 0 0 0 NA NA
## vxreg 0 0 0 0 NA NA
## skew 0 0 1 1 NA NA
## shape 0 0 1 1 NA NA
## ghlambda 0 0 0 0 NA NA
## xi 0 0 0 0 NA NA#fit the model to the data
modelfit=ugarchfit(spec=model,data=returns)
modelfit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(14,0)
## Mean Model : ARFIMA(2,0,2)
## Distribution : sstd
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000424 0.000115 3.694424 0.000220
## ar1 0.028714 0.333601 0.086074 0.931407
## ar2 0.379502 0.217792 1.742496 0.081422
## ma1 -0.106715 0.330267 -0.323117 0.746606
## ma2 -0.446759 0.238626 -1.872216 0.061177
## omega 0.000017 0.000003 6.637722 0.000000
## alpha1 0.000000 0.014376 0.000000 1.000000
## alpha2 0.124089 0.025264 4.911631 0.000001
## alpha3 0.090595 0.023149 3.913513 0.000091
## alpha4 0.092997 0.024990 3.721423 0.000198
## alpha5 0.104300 0.024761 4.212223 0.000025
## alpha6 0.061165 0.022886 2.672616 0.007526
## alpha7 0.068429 0.023645 2.894055 0.003803
## alpha8 0.096837 0.025670 3.772333 0.000162
## alpha9 0.052590 0.021073 2.495599 0.012574
## alpha10 0.072525 0.024225 2.993823 0.002755
## alpha11 0.072022 0.022884 3.147225 0.001648
## alpha12 0.022210 0.020987 1.058290 0.289923
## alpha13 0.036099 0.019535 1.847881 0.064620
## alpha14 0.023143 0.019559 1.183250 0.236710
## skew 0.870222 0.021124 41.195553 0.000000
## shape 8.630207 1.297672 6.650531 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000424 0.000129 3.29899 0.000970
## ar1 0.028714 0.216027 0.13292 0.894256
## ar2 0.379502 0.114688 3.30898 0.000936
## ma1 -0.106715 0.213577 -0.49966 0.617316
## ma2 -0.446759 0.126705 -3.52599 0.000422
## omega 0.000017 0.000002 7.76844 0.000000
## alpha1 0.000000 0.017545 0.00000 1.000000
## alpha2 0.124089 0.029129 4.25998 0.000020
## alpha3 0.090595 0.021727 4.16974 0.000030
## alpha4 0.092997 0.025845 3.59828 0.000320
## alpha5 0.104300 0.023882 4.36731 0.000013
## alpha6 0.061165 0.023226 2.63348 0.008451
## alpha7 0.068429 0.024684 2.77223 0.005567
## alpha8 0.096837 0.027525 3.51820 0.000434
## alpha9 0.052590 0.020928 2.51295 0.011973
## alpha10 0.072525 0.025149 2.88376 0.003930
## alpha11 0.072022 0.022984 3.13356 0.001727
## alpha12 0.022210 0.023142 0.95973 0.337192
## alpha13 0.036099 0.018060 1.99882 0.045628
## alpha14 0.023143 0.021759 1.06361 0.287507
## skew 0.870222 0.020989 41.46165 0.000000
## shape 8.630207 1.277840 6.75375 0.000000
##
## LogLikelihood : 10562.67
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.2854
## Bayes -6.2453
## Shibata -6.2855
## Hannan-Quinn -6.2711
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.738 1.874e-01
## Lag[2*(p+q)+(p+q)-1][11] 11.264 1.024e-12
## Lag[4*(p+q)+(p+q)-1][19] 17.604 3.533e-03
## d.o.f=4
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.2877 0.5917
## Lag[2*(p+q)+(p+q)-1][41] 13.3219 0.9431
## Lag[4*(p+q)+(p+q)-1][69] 23.4313 0.9680
## d.o.f=14
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[15] 1.900 0.500 2.000 0.1680
## ARCH Lag[17] 2.027 1.496 1.887 0.5408
## ARCH Lag[19] 2.843 2.483 1.802 0.6715
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: no.parameters>20 (not available)
## Individual Statistics:
## mu 0.71695
## ar1 0.08409
## ar2 0.04561
## ma1 0.09545
## ma2 0.04949
## omega 0.39837
## alpha1 0.38957
## alpha2 0.32477
## alpha3 0.22171
## alpha4 0.06361
## alpha5 0.30666
## alpha6 0.06821
## alpha7 0.19388
## alpha8 0.42707
## alpha9 0.06410
## alpha10 0.06916
## alpha11 0.23119
## alpha12 0.09882
## alpha13 0.13785
## alpha14 0.09918
## skew 0.31054
## shape 0.47752
##
## Asymptotic Critical Values (10% 5% 1%)
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.8962 5.803e-02 *
## Negative Sign Bias 0.6498 5.159e-01
## Positive Sign Bias 2.7377 6.220e-03 ***
## Joint Effect 30.4067 1.133e-06 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 49.09 0.0001782
## 2 30 56.61 0.0016004
## 3 40 65.99 0.0044346
## 4 50 69.67 0.0276760
##
##
## Elapsed time : 6.710096plot(modelfit, which = 10)
plot(modelfit, which = 11)
Looking at the ACF of squared standardized residuals, it seems that the
values stay within the bands, meaning dynamics have been wiped out. This
model is good for fitting the data.
#After comparing AIC/BICs for various GARCH combinations, I decided to go with an GARCH(2,1) model, which had the lowest value.
model1=ugarchspec(
variance.model = list(model = "sGARCH", garchOrder = c(2, 1)),
mean.model = list(armaOrder = c(2, 2), include.mean = TRUE),
distribution.model = "sstd")
#sanity check: explore the model parameters
model1@model$pars## Level Fixed Include Estimate LB UB
## mu 0 0 1 1 NA NA
## ar1 0 0 1 1 NA NA
## ar2 0 0 1 1 NA NA
## ma1 0 0 1 1 NA NA
## ma2 0 0 1 1 NA NA
## arfima 0 0 0 0 NA NA
## archm 0 0 0 0 NA NA
## mxreg 0 0 0 0 NA NA
## omega 0 0 1 1 NA NA
## alpha1 0 0 1 1 NA NA
## alpha2 0 0 1 1 NA NA
## beta1 0 0 1 1 NA NA
## gamma 0 0 0 0 NA NA
## eta1 0 0 0 0 NA NA
## eta2 0 0 0 0 NA NA
## delta 0 0 0 0 NA NA
## lambda 0 0 0 0 NA NA
## vxreg 0 0 0 0 NA NA
## skew 0 0 1 1 NA NA
## shape 0 0 1 1 NA NA
## ghlambda 0 0 0 0 NA NA
## xi 0 0 0 0 NA NA#fit the model to the data
modelfit1=ugarchfit(spec=model1,data=returns)
modelfit1##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(2,1)
## Mean Model : ARFIMA(2,0,2)
## Distribution : sstd
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000414 0.000101 4.082633 0.000045
## ar1 0.075122 0.347045 0.216460 0.828629
## ar2 0.364454 0.232917 1.564738 0.117644
## ma1 -0.155823 0.344871 -0.451829 0.651392
## ma2 -0.426434 0.257215 -1.657887 0.097340
## omega 0.000002 0.000006 0.293893 0.768839
## alpha1 0.000000 0.014750 0.000032 0.999974
## alpha2 0.114612 0.089699 1.277738 0.201342
## beta1 0.876121 0.087667 9.993787 0.000000
## skew 0.866983 0.012883 67.299125 0.000000
## shape 8.626200 2.659506 3.243535 0.001181
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000414 0.000879 0.470837 0.637757
## ar1 0.075122 0.156581 0.479760 0.631398
## ar2 0.364454 0.902869 0.403662 0.686461
## ma1 -0.155823 0.154911 -1.005884 0.314471
## ma2 -0.426434 0.902482 -0.472512 0.636561
## omega 0.000002 0.000094 0.017788 0.985808
## alpha1 0.000000 0.075751 0.000006 0.999995
## alpha2 0.114612 1.516198 0.075592 0.939744
## beta1 0.876121 1.457721 0.601021 0.547826
## skew 0.866983 0.276464 3.135972 0.001713
## shape 8.626200 48.751382 0.176943 0.859553
##
## LogLikelihood : 10563.23
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.2923
## Bayes -6.2723
## Shibata -6.2923
## Hannan-Quinn -6.2851
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.598 2.061e-01
## Lag[2*(p+q)+(p+q)-1][11] 11.520 1.323e-13
## Lag[4*(p+q)+(p+q)-1][19] 18.089 2.305e-03
## d.o.f=4
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.005113 0.9430
## Lag[2*(p+q)+(p+q)-1][8] 2.196073 0.8277
## Lag[4*(p+q)+(p+q)-1][14] 3.905436 0.8891
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 0.002388 0.500 2.000 0.9610
## ARCH Lag[6] 1.155457 1.461 1.711 0.7036
## ARCH Lag[8] 1.802242 2.368 1.583 0.7812
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 97.0513
## Individual Statistics:
## mu 0.73947
## ar1 0.09432
## ar2 0.05534
## ma1 0.10374
## ma2 0.06097
## omega 20.67564
## alpha1 0.20342
## alpha2 0.18104
## beta1 0.17380
## skew 0.38641
## shape 0.49342
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 2.49 2.75 3.27
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.9387 5.262e-02 *
## Negative Sign Bias 0.5138 6.074e-01
## Positive Sign Bias 2.7341 6.288e-03 ***
## Joint Effect 29.9240 1.432e-06 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 41.79 0.001892
## 2 30 54.37 0.002931
## 3 40 56.81 0.032544
## 4 50 74.62 0.010622
##
##
## Elapsed time : 1.991244plot(modelfit1, which = 10)
plot(modelfit1, which = 11)
The GARCH model fits equally well since there are generally no spikes.
It also has a slightly lower AIC/BIC. I will prefer the GARCH model due
to having less parameters compared to the ARCH model.
2) Problem 14.4
#forecast
modelfor = ugarchforecast(modelfit1, data = NULL, n.ahead = 2, n.roll = 0, out.sample = 0)
plot(modelfor, which = 1)
#get forecasted volatility
forecasted_volatility <- sigma(modelfor)
#get 95% prediction intervals for 1 and 2 step ahead forecasts
lower_1 <- -1.96 * forecasted_volatility[1]
upper_1 <- 1.96 * forecasted_volatility[1]
lower_2 <- -1.96 * forecasted_volatility[2]
upper_2 <- 1.96 * forecasted_volatility[2]
#print 95% prediction intervals
cat("One-Step Ahead - Lower Bound:", lower_1, "Upper Bound:", upper_1, "\n")## One-Step Ahead - Lower Bound: -0.01713364 Upper Bound: 0.01713364cat("Two-Step Ahead - Lower Bound:", lower_2, "Upper Bound:", upper_2, "\n")## Two-Step Ahead - Lower Bound: -0.01634723 Upper Bound: 0.016347233) Problem 14.5
cpi <- read_excel("Chapter14_exercises_data.xls", sheet = "Exercise 5a")
gdp <- read_excel("Chapter14_exercises_data.xls", sheet = "Exercise 5b")
inf <- cpi$`INFLATION RATE`[2:nrow(cpi)]
gdp_grow <- gdp$`GDP GROWTH`[2:nrow(gdp)]
mean(inf)## [1] 0.3009486mean(gdp_grow)## [1] 0.7842689plot(inf, type="l")
plot(gdp_grow, type = "l")
adf.test(cpi$`INFLATION RATE`[2:nrow(cpi)])##
## Augmented Dickey-Fuller Test
##
## data: cpi$`INFLATION RATE`[2:nrow(cpi)]
## Dickey-Fuller = -3.9489, Lag order = 9, p-value = 0.01155
## alternative hypothesis: stationaryadf.test(gdp$`GDP GROWTH`[2:nrow(gdp)])## Warning in adf.test(gdp$`GDP GROWTH`[2:nrow(gdp)]): p-value smaller than
## printed p-value##
## Augmented Dickey-Fuller Test
##
## data: gdp$`GDP GROWTH`[2:nrow(gdp)]
## Dickey-Fuller = -6.4572, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationaryauto.arima(cpi$"INFLATION RATE"[2:nrow(cpi)]) #ARIMA order (4,1,1)## Series: cpi$"INFLATION RATE"[2:nrow(cpi)]
## ARIMA(4,1,1)
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1
## 0.2872 -0.0455 0.0071 -0.0227 -0.8607
## s.e. 0.0715 0.0509 0.0495 0.0525 0.0626
##
## sigma^2 = 0.078: log likelihood = -111.72
## AIC=235.45 AICc=235.55 BIC=263.5auto.arima(gdp$"GDP GROWTH"[2:nrow(gdp)]) #ARIMA order (2,0,2)## Series: gdp$"GDP GROWTH"[2:nrow(gdp)]
## ARIMA(2,0,2) with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 ma2 mean
## 1.2589 -0.6531 -0.9388 0.4793 0.7814
## s.e. 0.1742 0.1593 0.2075 0.1654 0.0763
##
## sigma^2 = 0.8324: log likelihood = -348
## AIC=708 AICc=708.33 BIC=729.46diff_inf<- diff(inf)#Model for inflation rate
model <- ugarchspec(
variance.model = list(model = "sGARCH", garchOrder = c(3, 3)),
mean.model = list(armaOrder = c(4, 1), include.mean = TRUE),
distribution.model = "sstd"
)
# Sanity check: explore the model parameters
print(model@model$pars)## Level Fixed Include Estimate LB UB
## mu 0 0 1 1 NA NA
## ar1 0 0 1 1 NA NA
## ar2 0 0 1 1 NA NA
## ar3 0 0 1 1 NA NA
## ar4 0 0 1 1 NA NA
## ma1 0 0 1 1 NA NA
## arfima 0 0 0 0 NA NA
## archm 0 0 0 0 NA NA
## mxreg 0 0 0 0 NA NA
## omega 0 0 1 1 NA NA
## alpha1 0 0 1 1 NA NA
## alpha2 0 0 1 1 NA NA
## alpha3 0 0 1 1 NA NA
## beta1 0 0 1 1 NA NA
## beta2 0 0 1 1 NA NA
## beta3 0 0 1 1 NA NA
## gamma 0 0 0 0 NA NA
## eta1 0 0 0 0 NA NA
## eta2 0 0 0 0 NA NA
## delta 0 0 0 0 NA NA
## lambda 0 0 0 0 NA NA
## vxreg 0 0 0 0 NA NA
## skew 0 0 1 1 NA NA
## shape 0 0 1 1 NA NA
## ghlambda 0 0 0 0 NA NA
## xi 0 0 0 0 NA NA# Fit the model to the data
modelfit1 <- ugarchfit(spec = model, data = diff_inf)
print(modelfit1)##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(3,3)
## Mean Model : ARFIMA(4,0,1)
## Distribution : sstd
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001238 0.001359 0.911414 0.362077
## ar1 0.196110 0.060187 3.258320 0.001121
## ar2 -0.004106 0.049146 -0.083547 0.933417
## ar3 -0.042455 0.046840 -0.906382 0.364734
## ar4 -0.019913 0.043901 -0.453588 0.650125
## ma1 -0.838828 0.043947 -19.087256 0.000000
## omega 0.008984 0.011889 0.755648 0.449860
## alpha1 0.197500 0.072165 2.736772 0.006205
## alpha2 0.168374 0.344144 0.489255 0.624661
## alpha3 0.124121 0.233845 0.530782 0.595570
## beta1 0.000000 1.771586 0.000000 1.000000
## beta2 0.320327 0.786156 0.407460 0.683670
## beta3 0.120003 0.262842 0.456561 0.647987
## skew 1.086276 0.049806 21.809933 0.000000
## shape 4.698211 0.676905 6.940723 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001238 0.001565 0.791306 0.428766
## ar1 0.196110 0.102226 1.918398 0.055061
## ar2 -0.004106 0.077936 -0.052684 0.957984
## ar3 -0.042455 0.076323 -0.556255 0.578036
## ar4 -0.019913 0.082114 -0.242504 0.808390
## ma1 -0.838828 0.116500 -7.200259 0.000000
## omega 0.008984 0.051475 0.174529 0.861450
## alpha1 0.197500 0.103260 1.912642 0.055794
## alpha2 0.168374 1.496250 0.112531 0.910403
## alpha3 0.124121 1.053896 0.117773 0.906247
## beta1 0.000000 7.290230 0.000000 1.000000
## beta2 0.320327 3.509238 0.091281 0.927269
## beta3 0.120003 0.884758 0.135634 0.892111
## skew 1.086276 0.084951 12.787098 0.000000
## shape 4.698211 1.629176 2.883795 0.003929
##
## LogLikelihood : 51.65621
##
## Information Criteria
## ------------------------------------
##
## Akaike -0.0924495
## Bayes -0.0040039
## Shibata -0.0931475
## Hannan-Quinn -0.0584586
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 6.373e-05 9.936e-01
## Lag[2*(p+q)+(p+q)-1][14] 1.195e+01 9.639e-11
## Lag[4*(p+q)+(p+q)-1][24] 2.828e+01 3.063e-06
## d.o.f=5
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 11.36 0.0007519
## Lag[2*(p+q)+(p+q)-1][17] 17.38 0.0219348
## Lag[4*(p+q)+(p+q)-1][29] 19.65 0.1496115
## d.o.f=6
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[7] 0.08025 0.500 2.000 0.7770
## ARCH Lag[9] 0.84418 1.485 1.796 0.8114
## ARCH Lag[11] 0.90917 2.440 1.677 0.9505
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 4.932
## Individual Statistics:
## mu 0.04198
## ar1 0.24350
## ar2 0.98424
## ar3 0.21016
## ar4 0.10888
## ma1 0.35696
## omega 0.13734
## alpha1 0.19043
## alpha2 0.16278
## alpha3 0.33587
## beta1 0.21233
## beta2 0.30283
## beta3 0.28592
## skew 0.06654
## shape 0.26498
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 3.26 3.54 4.07
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.4630 0.143874
## Negative Sign Bias 0.4078 0.683530
## Positive Sign Bias 3.0559 0.002319 ***
## Joint Effect 9.5390 0.022921 **
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 19.79 0.4075
## 2 30 27.74 0.5316
## 3 40 37.62 0.5329
## 4 50 61.92 0.1018
##
##
## Elapsed time : 1.085289plot(modelfit1, which = 10)
plot(modelfit1, which = 11)
modelfor = ugarchforecast(modelfit1, data = NULL, n.ahead = 2, n.roll = 0, out.sample = 0)
plot(modelfor, which = 1)
#get forecasted volatility
forecasted_volatility <- sigma(modelfor)
#One step ahead volatility forecast
lower_1 <- -1.96 * forecasted_volatility[1]
upper_1 <- 1.96 * forecasted_volatility[1]
#print 95% prediction intervals
cat("One-Step Ahead - Lower Bound:", lower_1, "Upper Bound:", upper_1, "\n")## One-Step Ahead - Lower Bound: -0.7550893 Upper Bound: 0.7550893#Model for GDP Growth
model_2 <- ugarchspec(
variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(2, 2), include.mean = TRUE),
distribution.model = "sstd"
)
# Sanity check: explore the model parameters
print(model_2@model$pars)## Level Fixed Include Estimate LB UB
## mu 0 0 1 1 NA NA
## ar1 0 0 1 1 NA NA
## ar2 0 0 1 1 NA NA
## ma1 0 0 1 1 NA NA
## ma2 0 0 1 1 NA NA
## arfima 0 0 0 0 NA NA
## archm 0 0 0 0 NA NA
## mxreg 0 0 0 0 NA NA
## omega 0 0 1 1 NA NA
## alpha1 0 0 1 1 NA NA
## beta1 0 0 1 1 NA NA
## gamma 0 0 0 0 NA NA
## eta1 0 0 0 0 NA NA
## eta2 0 0 0 0 NA NA
## delta 0 0 0 0 NA NA
## lambda 0 0 0 0 NA NA
## vxreg 0 0 0 0 NA NA
## skew 0 0 1 1 NA NA
## shape 0 0 1 1 NA NA
## ghlambda 0 0 0 0 NA NA
## xi 0 0 0 0 NA NAmodelfit_2 <- ugarchfit(spec = model, data = gdp_grow)
print(modelfit_2)##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(3,3)
## Mean Model : ARFIMA(4,0,1)
## Distribution : sstd
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.768010 0.173963 4.41478 0.000010
## ar1 1.192082 0.229627 5.19139 0.000000
## ar2 -0.071943 0.160867 -0.44722 0.654714
## ar3 -0.342749 0.108333 -3.16386 0.001557
## ar4 0.134246 0.068059 1.97248 0.048554
## ma1 -0.832107 0.246534 -3.37523 0.000738
## omega 0.082043 0.215821 0.38014 0.703839
## alpha1 0.271551 0.118106 2.29921 0.021493
## alpha2 0.212252 0.362164 0.58607 0.557831
## alpha3 0.000000 0.672449 0.00000 1.000000
## beta1 0.000000 0.280433 0.00000 1.000000
## beta2 0.000000 0.657710 0.00000 1.000000
## beta3 0.451791 0.116947 3.86322 0.000112
## skew 1.018403 0.128628 7.91742 0.000000
## shape 7.818062 4.907829 1.59298 0.111165
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.768010 1.02994 0.745683 0.45586
## ar1 1.192082 0.14201 8.394126 0.00000
## ar2 -0.071943 0.69750 -0.103145 0.91785
## ar3 -0.342749 0.23745 -1.443454 0.14889
## ar4 0.134246 0.18383 0.730264 0.46523
## ma1 -0.832107 0.68262 -1.218988 0.22285
## omega 0.082043 1.35318 0.060629 0.95165
## alpha1 0.271551 0.26193 1.036711 0.29987
## alpha2 0.212252 2.32039 0.091473 0.92712
## alpha3 0.000000 4.35927 0.000000 1.00000
## beta1 0.000000 1.90540 0.000000 1.00000
## beta2 0.000000 4.61083 0.000000 1.00000
## beta3 0.451791 0.28323 1.595127 0.11068
## skew 1.018403 0.60007 1.697139 0.08967
## shape 7.818062 20.95117 0.373156 0.70903
##
## LogLikelihood : -319.1496
##
## Information Criteria
## ------------------------------------
##
## Akaike 2.5314
## Bayes 2.7346
## Shibata 2.5254
## Hannan-Quinn 2.6131
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.5389 0.4629
## Lag[2*(p+q)+(p+q)-1][14] 5.6791 0.9995
## Lag[4*(p+q)+(p+q)-1][24] 11.9004 0.5632
## d.o.f=5
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.0158 0.9000
## Lag[2*(p+q)+(p+q)-1][17] 6.1031 0.7881
## Lag[4*(p+q)+(p+q)-1][29] 11.1602 0.7980
## d.o.f=6
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[7] 1.202 0.500 2.000 0.2730
## ARCH Lag[9] 2.467 1.485 1.796 0.4268
## ARCH Lag[11] 3.157 2.440 1.677 0.5671
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 3.0389
## Individual Statistics:
## mu 0.19668
## ar1 0.25910
## ar2 0.40438
## ar3 0.42239
## ar4 0.43056
## ma1 0.15904
## omega 0.52367
## alpha1 0.29816
## alpha2 0.38881
## alpha3 0.25425
## beta1 0.61003
## beta2 0.50725
## beta3 0.69710
## skew 0.03317
## shape 0.26153
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 3.26 3.54 4.07
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.90766 0.3649
## Negative Sign Bias 0.03468 0.9724
## Positive Sign Bias 0.59490 0.5524
## Joint Effect 3.39834 0.3342
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 11.00 0.9238
## 2 30 21.45 0.8421
## 3 40 44.18 0.2619
## 4 50 37.89 0.8752
##
##
## Elapsed time : 0.4812002plot(modelfit_2, which = 10)
plot(modelfit_2, which = 11)
modelfor = ugarchforecast(modelfit1, data = NULL, n.ahead = 2, n.roll = 0, out.sample = 0)
plot(modelfor, which = 1)
forecasted_volatility <- sigma(modelfor)
#One step ahead volatility forecast
lower_1 <- -1.96 * forecasted_volatility[1]
upper_1 <- 1.96 * forecasted_volatility[1]
#95% prediction intervals
cat("One-Step Ahead - Lower Bound:", lower_1, "Upper Bound:", upper_1, "\n")## One-Step Ahead - Lower Bound: -0.7550893 Upper Bound: 0.75508934. Problem 12.2
(a)
gasoline_ts <- ts(us_gasoline$Barrels, frequency=52, start = c(1991, 6))# lambda of Box-Cox transformation
lambda_gasoline <- BoxCox.lambda(gasoline_ts)
# select the number of Fourier pairs
min.AIC <- Inf
K_min.Aic <- 0
for(num in c(1:6)){
gasoline_tslm <- tslm(
gasoline_ts ~ trend + fourier(gasoline_ts, K = num),
lambda = lambda_gasoline
)
AIC <- CV(gasoline_tslm)["AIC"]
if(AIC < min.AIC){
min.AIC <- AIC
K_min.Aic <- num
}
}
# Make harmonic regression model
gasoline_tslm <- tslm(
gasoline_ts ~ trend + fourier(gasoline_ts, K = K_min.Aic),
lambda = lambda_gasoline
)
autoplot(gasoline_ts) + autolayer(gasoline_tslm$fitted.values)
The dynamic harmonic regression model fits the data better than the harmonic regression model because it is able to capture changes in the seasonality and trend over time.
(b)
# checkresiduals(gasoline_autoarima)
checkresiduals(gasoline_tslm)
##
## Breusch-Godfrey test for serial correlation of order up to 104
##
## data: Residuals from Linear regression model
## LM test = 932.9, df = 104, p-value < 2.2e-16The residuals from the dynamic harmonic regression are almost all positive and there are several significant lags in the ACF plot of the residuals, but the autocorrelation is generally quite small. The residuals do not have constant volatility over time. For the harmonic regression model, the residuals are clearly not white noise as they display both trend and changing volatility. The ACFs decay extremely slowly starting from around 0.7.
(c)
We could also try models like ARIMA and ETS, although they might not perform as well on high frequency like weekly data with complex seasonal patterns. We generally used them for monthly or quarterly data.
# Fit ARIMA model
fit_arima <- auto.arima(gasoline_ts)
summary(fit_arima)## Series: gasoline_ts
## ARIMA(0,1,1)(1,0,0)[52]
##
## Coefficients:
## ma1 sar1
## -0.7440 0.2618
## s.e. 0.0206 0.0294
##
## sigma^2 = 0.07239: log likelihood = -144.86
## AIC=295.72 AICc=295.74 BIC=311.35
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.003688876 0.268748 0.2067944 -0.04355785 2.480655 0.7176306
## ACF1
## Training set -0.01839917checkresiduals(fit_arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(1,0,0)[52]
## Q* = 330.24, df = 102, p-value < 2.2e-16
##
## Model df: 2. Total lags used: 104# Fit ETS model
fit_ets <- ets(gasoline_ts)## Warning in ets(gasoline_ts): I can't handle data with frequency greater than
## 24. Seasonality will be ignored. Try stlf() if you need seasonal forecasts.summary(fit_ets)## ETS(A,N,N)
##
## Call:
## ets(y = gasoline_ts)
##
## Smoothing parameters:
## alpha = 0.3251
##
## Initial states:
## l = 6.7151
##
## sigma: 0.2772
##
## AIC AICc BIC
## 6299.071 6299.089 6314.705
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.003778683 0.2770348 0.2132618 -0.04538429 2.560032 0.7400743
## ACF1
## Training set -0.002521048checkresiduals(fit_ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(A,N,N)
## Q* = 455.07, df = 104, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 104autoplot(gasoline_ts, ylab = "Gas Supply (Weekly)",main= "Other Models") + autolayer(fitted(fit_arima)) + autolayer(fitted(fit_ets))
5) 12.3
# experiment with retail data.
# create retail data
retail <- read.xlsx("retail.xlsx",
sheetIndex = 1,
startRow = 2)
retail.ts <- ts(retail[,"A3349873A"],
frequency=12,
start=c(1982,4))
# fit model
retail_nnetar <- nnetar(retail.ts)
checkresiduals(retail_nnetar)
##
## Ljung-Box test
##
## data: Residuals from NNAR(3,1,2)[12]
## Q* = 311.51, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24autoplot(retail.ts) + autolayer(retail_nnetar$fitted)
# Produce forecasts
fc_retail_nnetar <- forecast(retail_nnetar, h = 36)
autoplot(fc_retail_nnetar)
# test accuracy using future data.
retail.new <- read.xlsx("8501011.xlsx",
sheetName = "Data1",
startRow = 10)
retail.new.ts <- ts(retail.new[, "A3349873A"],
start = c(1982, 4),
frequency = 12)
retail.new.test <- subset(
retail.new.ts,
start = length(retail.ts) + 1
)
# test accuracy using future data.
accuracy(fc_retail_nnetar, retail.new.test)## ME RMSE MAE MPE MAPE MASE
## Training set 0.01200432 19.38566 13.77693 -0.7347354 6.078575 0.7276128
## Test set -51.12171115 69.70154 64.71798 -12.0187699 14.337653 3.4180073
## ACF1 Theil's U
## Training set 0.5539058 NA
## Test set 0.6514395 1.434885Even the the residuals aren’t like white nice and were slightly right skewed, it looks like neural network model fitted well to the data.
# experiment with ibmclose data.
ibmclose_nnetar <- nnetar(ibmclose)
# Produce forecasts
autoplot(ibmclose) + autolayer(ibmclose_nnetar$fitted)
checkresiduals(ibmclose_nnetar)
##
## Ljung-Box test
##
## data: Residuals from NNAR(1,1)
## Q* = 15.703, df = 10, p-value = 0.1085
##
## Model df: 0. Total lags used: 10fc_ibmclose_nnetar <- forecast(ibmclose_nnetar, h=36)
autoplot(fc_ibmclose_nnetar)
The fitted value looks pretty good and residuals looks like white noise. However, even neural network method yielded naive-method like result in forecast. It looked like there wasn’t any rule in lagged values.
# experiment with usmelec data.
usmelec_nnetar <- nnetar(usmelec)
checkresiduals(usmelec_nnetar)
##
## Ljung-Box test
##
## data: Residuals from NNAR(3,1,2)[12]
## Q* = 215.96, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24autoplot(usmelec) + autolayer(usmelec_nnetar$fitted)
fc_usmelec_nnetar <- forecast(
usmelec_nnetar, h = 12*4
)
autoplot(fc_usmelec_nnetar)
It looked like neural network model was fitted well to the data.